Computing apparatus and computing method

ABSTRACT

An object of the invention is to provide a computing technology which can operate at room temperature and have a sufficient performance for combinatorial optimization problems that need an exhaustive search. In a local-field response method in which spins being variables respond to local effective magnetic fields, a time axis is discretely treated. When spins respond to effective magnetic fields, the effective magnetic fields are determined sequentially from the site having the small magnitude of a spin, and spins respond to the fields in order. When the sign of a spin is inverted, the information is reflected in the subsequent process of determining the effective magnetic fields for other sites. Thus, a many-body effect due to quantum entanglement is phenomenologically incorporated.

CROSS-REFERENCE TO RELATED APPLICATION

The present application claims priority from Japanese application JP2017-204990, filed on Oct. 24, 2017, the content of which is herebyincorporated by reference into this application.

BACKGROUND OF THE INVENTION 1. Field of the Invention

The present invention relates to a computing technology which is able toperform computation at a high speed with respect to combinatorialoptimization problems that need an exhaustive search.

2. Related Art

As being representative as the word “IoT” (Internet of Things), variousthings are connected to the Internet in the present days; information iscollected from the things; and the things are controlled by thecollected information. In the control, an optimal solution is found fromamong many choices and is performed. Extremely speaking, the informationtechnology in the present days can be said to search an optimalsolution.

In this situation, a quantum annealing, or adiabatic quantum computing,has recently received attention. In this method, a problem is set suchthat the ground state of a certain physical system is a solution and thesolution is obtained through finding the ground state. Let theHamiltonian of a physical system in which a problem is set be Ĥ_(p). Atthe beginning of computation, however, Hamiltonian is not Ĥ_(p) but Ĥ₀for which the ground state is easily prepared. Next, the Hamiltonian istransformed from Ĥ₀ to Ĥ_(p) with a sufficient period of time. When thetransformation takes a sufficient period of time, the system continuesto stay in the ground state, and the ground state (solution state) forHamiltonian Ĥ_(p) is finally obtained. This is a principle of thequantum annealing.

A ground state-searching method in which a physical system called Isingspin glass is used can be applicable even to a problem called NP-hard.Meanwhile, a highly difficult problem in combinatorial optimizationproblems belongs to the NP-hard. Moreover, problems classified into “P”and problems classified into “NP” in the computational complexity theoryall can be transformed to an NP-hard problem. Therefore, if the quantumannealing is applied to the Ising spin glass system, almost all thecombinatorial optimization problems can be solved, and the mostimportant issue of the information technology is solved.

Another reason why the quantum annealing receives attention isrobustness with respect to decoherence. In a quantum computer, quantumcoherence must be preserved over a computing time. However, in thequantum annealing, the condition is relaxed. If the ground state ismaintained, a correct solution is obtained. The quantum coherence is notnecessarily to be maintained. In a current technology level, it isdifficult to establish a pure quantum system. Therefore, the quantumcoherence is hardly kept over the computing time. For this reason, thequantum annealing has received attention. However, there is a drawbackeven in the quantum annealing. The quantum annealing can only berealized in a superconducting magnetic flux qubit system as it stands,and thus, a cryogenic cooling apparatus is needed. The need for anextremely low temperature is an issue of achieving a practical computer.

To solve the issue, the local-field response method was proposed asdescribed in the following (PTLs 1-3 and NPL 1). First, let us reviewthe quantum annealing again. The concept of the annealing liesregardless of quantum or classical by its nature. The quantum annealingis devised to improve the performance of the classic annealing usingquantum properties. The reason why the quantum annealing does notrequire the quantum coherence over the computing time and only requiresthe ground state being maintained comes from the wide concept of theannealing. As just described, the concept of the annealing is so widethat another method to use quantum properties may exist, different fromthe quantum annealing.

The local-field response method has been disclosed from such a point ofview. In this method, similarly to the quantum annealing, a transversefield is applied to a spin system at time t=t₀, which is a computingdevice, and the magnetic field is gradually decreased to obtain asolution at time t=τ. The computing device itself is classical, andquantum-mechanical information is added to the response of spins to themagnetic field. The method works on a classical machine at roomtemperature, and therefore, it can solve the issue of cryogenic coolingin the quantum annealing. In PTLs 1-3 and NPL 1, quantum effects wereintroduced in the response function, where the response function whichhas an average quantum effect was determined empirically or based on theresults solved for similar problems. In NPL 2, solution accuracy wasimproved compared to the methods in PTLs 1-3 and NPL 1 byphenomenologically incorporating the properties of the linearsuperposition in quantum mechanics. However, the property of quantumentanglement which is another important property in quantum mechanicshas not been sufficiently incorporated.

CITATION LIST Patent Literature

-   PTL 1: WO 2015/118639-   PTL 2: WO 2016/157333-   PTL 3: WO 2016/194221

Non-Patent Literature

-   NPL 1: T. Tomaru, “Quasi-Adiabatic Quantum Computing Treated with    c-Numbers Using the Local-Field Response,” J. Phys. Soc. Jpn. 85,    034802 (2016).-   NPL 2: T. Tomaru, “Improved Local-field response Method Working as    Quasi-Quantum Annealing,” J. Phys. Soc. Jpn. 86, 054801 (2017).

SUMMARY OF INVENTION Technical Problem

As described above, the quantum annealing requires a cryogenic coolingapparatus to use a superconducting magnetic flux qubit system.Meanwhile, the local-field response method operates at room temperature,but quantum entanglement which is an important property in quantummechanics is not sufficiently incorporated, which limits theperformance. Thus, an object of the invention is to provide a computingapparatus and a computing program which can operate at room temperatureand have a sufficient performance for a difficult assignment such as anissue that needs an exhaustive search.

Solution to Problem

In a local-field response method in which spins being variables respondto local effective magnetic fields, a time axis is discretely treated;When spins respond to effective magnetic fields, the effective magneticfields are determined sequentially from the site having the smallmagnitude of a spin, and the spins respond to the fields in order. Whenthe sign of the spin is inverted, the information is reflected in thesubsequent process of determining the effective magnetic fields forother sites. Thus, a many-body effect due to quantum entanglement isphenomenologically incorporated. More specific descriptions are providedin the following.

A specific figure is a computing apparatus which includes a computingunit, a storage unit, and a control unit, and performs computation underthe control of the control unit while transferring data between thestorage unit and the computing unit, or a computing method using thecomputing apparatus. N variables s_(j) ^(z) (j=1, 2, . . . , N) take arange of −1≤s_(j) ^(z)≤1, and an assignment is set with coefficientsg_(j) indicating local terms and coefficients J_(kj) (k, j=1, 2, . . . ,N) indicating inter-variable interactions. Time is divided into m, andthe computing unit discretely performs computation from t=t₀ (t₀=0) tot_(m) (t_(m)≤τ). Herein, N and m are natural numbers.

Variables B_(eff,j) ^(z)(t_(i)) and s_(j) ^(z)(t_(i)) at each time t_(i)(i=1, 2, . . . , m) are determined in this order. B_(eff,j) ^(z)(t_(i))is a function of s_(k) ^(z)(t_(i−1)), J_(kj), g_(j), and t_(i). s_(j)^(z)(t_(i)) is a function of B_(eff,j) ^(z)(t_(i)) and t_(i). Initialvalues at time t₀ are set as B_(j) ^(z)(t₀)=0 and s_(j) ^(z)(t₀)=0. Fordetermining B_(eff,j) ^(z)(t_(i)) and s_(j) ^(z)(t_(i)) at time t_(i)(i=1, 2, . . . , m), first, s_(j) ^(z)(t_(i−1)) are put in descendingorder such that |s_(m1) ^(z)(t_(i−1))|≤|s_(m2) ^(z)(t_(i−1))|≤s_(m3)^(z)(t_(i−1))|≤ . . . ≤|s_(mN) ^(z)(t_(i−1))|. Then, B_(eff,m1)^(z)(t_(i)) and s_(m1) ^(z)(t_(i)) at site m₁ are determined at thefirst time, and s_(m1) ^(z)(t_(i−1)) is set to be s_(m1)^(z)(t_(i−1))=sgn(s_(m1) ^(z)(t_(i)))|s_(m1) ^(z)(t_(i−1))|. Next,B_(eff,m2) ^(z)(t_(i)) and s_(m2) ^(z)(t_(i)) at site m₂ are determinedand s_(m2) ^(z)(t_(i−1)) is set to be s_(m2) ^(z)(t_(i−1))=sgn(s_(m2)^(z)(t_(i)))|s_(m2) ^(z)(t_(i−1))|. The computation at site m₃ isperformed similarly, and the computation up to site m_(N) is performedsimilarly for the computation at time t_(i). As the time step progressesfrom t=t₀ to t=t_(m), the variable s_(j) ^(z) approaches −1 or 1, and asolution is determined as s_(j) ^(zfd)=−1 if s_(j) ^(z)<0 and as s_(j)^(zfd)=1 if s_(j) ^(z)>0.

As described, the same computation has been performed up to the sitem_(N) for the computation at time t_(i), but there is an option that thesame computation is performed up to m_(x) (herein, 0<x≤N), and thatafter site m_(x) all remaining sites are processed in an individual andparallel manner similarly to the original local-field response methodfor the computation at time t_(i). Here, x is a natural number.

Advantageous Effects of Invention

The present method operates on a classical machine. Extremely lowtemperature is not needed. In addition, there is no need to take quantumcoherence into consideration. As a result, usable resources areexpanded, and electric circuits can also be used. Moreover, solutionaccuracy is improved and a computation time is shortened throughphenomenologically incorporating the effect of quantum entanglement.Thanks to these properties, it is possible to realize a practicalcomputing apparatus which can solve difficult problems with highsolution accuracy.

BRIEF DESCRIPTION OF DRAWINGS

FIG. 1 is a diagram schematically illustrating a principle of anembodiment;

FIG. 2 is a flowchart illustrating an example of an algorithm inaccordance with the first embodiment;

FIG. 3 is a graph plotting specific values of the response functionr_(b);

FIG. 4 is a flowchart illustrating an example of an algorithm inaccordance with the second embodiment;

FIG. 5 is a flowchart illustrating an example of an algorithm inaccordance with the third embodiment;

FIG. 6 is a flowchart illustrating another example of the algorithm inaccordance with the third embodiment;

FIG. 7 is a flowchart illustrating an example of an algorithm inaccordance with the fourth embodiment;

FIG. 8 is a flowchart illustrating an example of an algorithm inaccordance with the fifth embodiment;

FIG. 9 is a diagram illustrating an example of an algorithm inaccordance with the sixth embodiment; and

FIG. 10 is a block diagram illustrating an example of a configuration ofa computing apparatus in accordance with the seventh embodiment.

DESCRIPTION OF EMBODIMENTS

Embodiments will be described in detail using the drawings. However, thecontent of the embodiments described below is not interpreted in a wayof limiting the invention. A person skilled in the art can easilyunderstand that the specific configuration may vary in a scope notdeparting from the idea and the spirit of the invention.

Portions having the same or similar functions in the configuration ofthe invention described below will be represented with the same symboland commonly used in different drawings, and the redundant descriptionwill be omitted.

In a case where there are elements having similar or correspondingfunctions, the elements may be attached to the same symbol withdifferent suffixes. However, in a case where there is no need todistinguish plural elements, the suffixes may be omitted.

First Embodiment

In the first embodiment, we will give the basic principle, starting froma quantum-mechanical description and transforming it into a classicalform.

FIG. 1 schematically show the principle of the embodiment. A basicframework is the same as the local-field response method disclosed inPTL 1 and NPL 1. A transverse field is applied at t=0 to make spinsdirected in one direction. Thereafter, the transverse field is graduallydecreased, and Hamiltonian is set to the problem Hamiltonian at t=τ.Spins time-evolve in response to the local effective magnetic fieldwhich is applied to each spin at each time.

Let the problem Hamiltonian and the Hamiltonian at t=0 be Eqs. (1) and(2), respectively.

$\begin{matrix}{{\hat{H}}_{p} = {{- {\sum\limits_{i > j}{J_{ij}{\hat{\sigma}}_{i}^{z}{\hat{\sigma}}_{j}^{z}}}} - {\sum\limits_{j}{g_{j}{\hat{\sigma}}_{j}^{z}}}}} & (1) \\{{\hat{H}}_{0} = {{- \gamma}{\sum\limits_{j}{\hat{\sigma}}_{j}^{z}}}} & (2)\end{matrix}$

Let the Hamiltonian at time t be Eq. (3).

$\begin{matrix}{{\hat{H}(t)} = {{\left( {1 - \frac{t}{\tau}} \right){\hat{H}}_{0}} + {\frac{t}{\tau}{\hat{H}}_{p}}}} & (3)\end{matrix}$

Herein, τ is a computation time. From an analogy with a one-spin system,the effective magnetic field at site j is given by B̂_(eff,j)=−∂Ĥ/∂σ̂_(j).

$\begin{matrix}{{{\hat{B}}_{{eff},j}(t)} = \left( {{\left( {1 - \frac{t}{\tau}} \right)\gamma},0,{\frac{t}{\tau}\left( {{\sum\limits_{k{({\neq j})}}{J_{kj}{\hat{\sigma}}_{k}^{z}}} + g_{j}} \right)}} \right)} & (4)\end{matrix}$

The local-field response method of this embodiment operates on aclassical machine in which expectation value <σ̂_(j)> is regarded as aspin variable. As seen in Eqs. (1) and (2), <σ̂_(j)> and <B̂_(eff,j)>consists of only x and z components. Therefore, a response functionr_(b)(t) is defined only by the x and z components as described by Eq.(5).

{circumflex over (σ)}_(j) ^(z)(t)

/

{circumflex over (σ)}_(j) ^(x)(t)

=r _(b)(t)·

{circumflex over (B)} _(eff,j) ^(z)(t)

/

{circumflex over (B)} _(eff,j) ^(x)(t)

  (5)

A spin direction is determined based on the response function. When thespin system is classical, the response of each spin is determined onlyby the effective magnetic field at each site, and the response functionis r_(b)(t)=1.

However, there is a non-local correlation (quantum entanglement) inquantum mechanics, and generally r_(b)(t)≠1. As mentioned already, Eq.(5) has been transformed to a classic form by taking the expectationvalue, but a quantum effect is incorporated in r_(b)(t)≠1. A value ofr_(b)(t) is determined empirically or through quantum-mechanical priorcalculations for similar problems. The quantum effect herein is anaverage one because it is empirically determined or it is based on theresults for similar problems. The local-field response method canincorporate quantum effects through various methods. The method throughr_(b)(t)≠1 is an example. Note that setting r_(b)(t)=1 without thequantum effect is also allowed in the local-field response method. Thelocal-field response method itself can operate even when the quantumeffect is not incorporated.

Four variables <σ̂_(j) ^(z)(t_(i))>, <σ̂_(j) ^(x)(t_(i))>, <B̂_(eff,j)^(z)(t_(i))>, and <B̂_(eff,j) ^(x)(t_(i))> in Eq. (5) are expectationvalues, and therefore, classical quantities. For this reason, let uschange the quantum-mechanical description to a classical-physics one,i.e., <σ̂_(j) ^(x)(t_(i))>→s_(j) ^(x)(t_(i)), <σ̂_(j) ^(z)(t_(i))>→s_(j)^(z)(t_(i)), <B̂_(eff,j) ^(x)(t_(i))>→B_(eff,j) ^(x)(t_(i)), <B̂_(eff,j)^(z)(t_(i))>→B_(eff,j) ^(z)(t_(i)). With the description change, Eq. (5)is modified to Eq. (6).

(s _(j) ^(z)(t)/s _(j) ^(x)(t))=r _(b)(t)·(B _(eff,j) ^(z)(t)/B _(eff,j)^(x)(t))  (6)

The time-evolution in the local-field response method is discrete.B_(eff,j) ^(z)(t_(i)) at time t_(i) is determined from s_(k)^(z)(t_(i−1)) at time t_(i−1) in accordance with Eq. (4). s_(j)^(z)(t_(i)) at time t_(i) is determined from B_(eff,j) ^(z)(t_(i)) attime t_(i) in accordance with Eq. (6). This procedure is repeated. FIG.2 summarizes the procedure as a flowchart.

Step 101 in FIG. 2 indicates a start point of the algorithm and setsinitial values. Step 102 a determines B_(eff,j) ^(z)(t_(i)) using s_(k)^(z)(t_(i−1)) at time t_(i−1), and determines the transverse fieldintensity B_(eff,j) ^(x)(t_(i)) depending on time. Step 103 determinestan θ=r_(b)(t)·B_(eff,j) ^(z)(t_(i))/B_(eff,j) ^(x)(t_(i)) correspondingto the spin direction using B_(eff,j) ^(z)(t_(i))/B_(eff,j) ^(x)(t_(i))and the response function r_(b)(t), and it determines s_(j)^(z)(t_(i))=r_(s)(t)·sin θ using a parameter r_(s)(t) representing themagnitude of the spin. r_(s)(t) is defined by r_(s)(t)²=s_(j)^(x)(t_(i))²=s_(j) ^(z)(t_(i))², and it is determined empirically orthrough prior calculations similarly to r_(b)(t). A specific examplewill be described in the following second paragraph. Steps 102 a and 103in FIG. 2 are a set of repeated computation. The set is repeated untilt=t_(m) (≤τ). At t=t_(m), a solution s_(j) ^(zfd) is determined as s_(j)^(zfd)=1 if s_(j) ^(z)(t)>0 and s_(j) ^(zfd)=−1 if s_(j) ^(z)(t_(m))<0(Steps 201 and 202). Herein, the reason why t_(m)≤τ is assumed is thatmany cases obtain converged solutions before t=τ in the time-evolution.

Step 103 can be generally written using a function f such as s_(j)^(z)(t_(i))=f(B_(eff,j) ^(z)(t_(i)), t_(i)), and f(B_(eff,j)^(z)(t_(i)), t_(i))=r_(s)(t)·sin{arctan(r_(b)(t)·B_(eff,j)^(z)(t_(k))/B_(eff,j) ^(x)(t_(k)))}. r_(s)(t) is a parameter expressingthe magnitude of the spin and satisfies 0≤r_(s)(t)≤1. Thus, −1≤s_(j)^(z)(t_(i))≤1. Here, when r_(b)(t)=1 and r_(s)(t)=1, then the system ispurely classical.

FIG. 3 shows an example of prior computations for the response functionr_(b)(t). The example shows a case where J_(ij) and g_(j) are determinedby uniform random numbers of [−5, 5] in an 8-bit system. The figureshows the results of 100 problems and contains 800 points (100problems×8 bits). Here, B_(zx)≡<B̂_(eff,j) ^(z)>/<B̂_(eff,j) ^(x)> ands_(zx)≡<σ̂_(j) ^(z)>/<σ̂_(j) ^(x)>. The points indicate values computedexactly quantum-mechanically. The response functions largely scatterowing to the non-local correlation in quantum mechanics. Open circlesindicate average values, where the horizontal axis is divided into 40regions and the points in each region are averaged. The averagedresponse function has smooth dependence on B_(zx). Because being smooth,the response function can be expressed with several parameters. NPL 1discloses a technique of expressing the smooth response function usingfour parameters. A solid line r_(b) ⁰(t) in FIG. 3 is obtained with thetechnique. Another parameter r_(s)(t) is also obtained from the samefour parameters.

Hitherto, an example of the response function has been described in FIG.3, and a method of determining r_(b) ⁰(t) and r_(s)(t) has beendescribed with reference to NPL 1. However, the method of determiningthe response functions r_(b) ⁰(t) and r_(s)(t) is not limited to theabove description, and various methods may be employed and the responsefunctions can also be empirically determined. In the followingembodiment, r_(b) ⁰(t) is not limited to the solid line in FIG. 3, butassumed to be a response function in which quantum effects areincorporated on average.

Here, similarly to r_(b)(t), r_(s)(t) is also possible to setr_(s)(t)=1. While the range of r_(b)(t) is −∞<r_(b)(t)<∞, r_(s)(t) is0≤r_(s)(t)≤1. Therefore, the influence of assuming r_(s)(t)=1 on a finalsolution is smaller than that of assuming r_(b)(t)=1. Thus, it iseffective to set r_(s)(t)=1 when prior information for determiningr_(s)(t) is insufficient.

Second Embodiment

The first embodiment has been described that quantum effects can beaveragely incorporated through r_(b)(t)≠1. However, quantum effectsdepend on problems and vary with time. Quantum effects cannot besufficiently incorporated only through averaged quantities. Thisembodiment describes a method of phenomenologically incorporating aquantum entanglement related-quantum effect in the formulation dependingon the spin state at each time.

The influence of quantum entanglement appears as a many-body effect.When quantum entanglement is large, if a certain spin is inverted (itssign is inverted), another spin is simultaneously inverted with highprobability. In the algorithm of FIG. 2, the effective magnetic fieldB_(eff,j) ^(z)(t_(i)) at t=t_(i) is calculated site by site using s_(j)^(z)(t_(i−1)) at t=t_(i−1). The calculation is performed independentlysite by site and it is in a one-body approximation. Therefore,simultaneous inversion of spins has not been sufficiently taken intoconsideration. For this reason, we take simultaneous inversion of spinsinto consideration in accordance with the algorithm in FIG. 4.

In the vicinity of the spin inversion, s_(j) ^(z)≈0. Therefore, the spinthat has a smaller value of |s_(j) ^(z)| has higher probability ofinversion. Thus, let us first put |s_(j) ^(z)(t_(i−1))| at each timet_(i−1) in descending order.

Let m₁, m₂, m₃, . . . be the site number of |s_(j) ^(z)(t_(i−1))| in thedescendent order as illustrated in Step 111 in FIG. 4. B_(eff,m1)^(z)(t_(i)) at site m₁ is calculated using Eq. (4) (Step 102 a), ands_(m1) ^(z)(t_(i)) is calculated using Eq. (6) (Step 103). If there isno sign inversion in the time evolution from s_(m1) ^(z)(t_(i−1)) tos_(m1) ^(z)(t_(i)), the next procedure is to compute B_(eff,m2)^(z)(t_(i)) and s_(m2) ^(z)(t_(i)) at site m₂. If the sign is inverted,s_(m1) ^(z)(t_(i−1)) is changed to −s_(m1) ^(z)(t_(i−1)), and afterthat, B_(eff,m2) ^(z)(t_(i)) and s_(m2) ^(z)(t_(i)) are computed. Inother words, the sign of s_(m1) ^(z)(t_(i−1)) is changed into the signof s_(m1) ^(z)(t_(i)) as s_(m1) ^(z)(t_(i−1))→sgn(s_(m1)^(z)(t_(i)))·|s_(m1) ^(z)(t₁₋₁)| (Step 112). Here, sgn denotes a signfunction which returns 1 or −1 in accordance with the sign of theargument.

According to this procedure, the change for site m₁ at time t_(i) isimmediately reflected on the time evolution for site m₂. When sites m₁and m₂ are entangled, if the spin at site m₁ is inverted, the spin atsite m₂ is inverted with high probability. Such an effect isincorporated as the change of the sign described herein. In other words,the effect of quantum entanglement is incorporated phenomenologically.

The computed result of s_(m2) ^(z)(t_(i)) is also processed similarly tos_(m1) ^(z)(t_(i)) i.e., s_(m2) ^(z)(t_(i−1))→sgn(s_(m2)^(z)(t_(i)))·|s_(m2) ^(z)(t_(i−1))| (Step 112). By repeating the similarprocess for site m₃ and the following sites, we obtain N components ofs_(j) ^(z)(t_(i)) (Step 113), and the process at time t_(i) iscompleted.

Once the process at time t_(i) is completed, the next procedure is theprocess at time t_(i+1), and the similar processes are repeated untiltime t=t_(m)≤τ to obtain the final solution. The final solution is s_(j)^(zfd)=1 if s_(j) ^(z)(t_(m))>0 and s_(j) ^(zfd)=−1 if s_(j)^(z)(t_(m))<0 similarly to the case of the first embodiment (Steps 201and 202).

The spin inversion is connected to a tunneling phenomenon in a viewpointof quantum mechanics. According to this interpretation, the treatment inthis embodiment is said to take multiple tunneling into consideration.

The energy value for the obtained solution is given by Eq. (7).

$\begin{matrix}{{H_{p}\left( t_{i} \right)} = {{- {\sum\limits_{k > j}{J_{kj}{s_{k}^{zfd}\left( t_{i} \right)}{s_{j}^{zfd}\left( t_{i} \right)}}}} - {\sum\limits_{j}{g_{j}{s_{j}^{zfd}\left( t_{i} \right)}}}}} & (7)\end{matrix}$

The local-field response method operates similarly to quantum annealing.The spin system in which the ground state is prepared at t=0time-evolves and comes to the ground state of the problem-embeddedsystem at t=τ ideally. The convergence of solutions is high according tothe method of this embodiment. The lowest energy state during thecomputation is the state at t=t_(m) with high probability (the groundstate with high probability). However, in the method of the firstembodiment, the convergence of the solution is low, and the lowestenergy state might be the state at t<t_(m). When the first embodiment isused, the energy needs to be calculated at each time using Eq. (7) andthe state corresponding to the lowest energy needs to be selected as thefinal solution. The amount of computation for that purpose is O(N²) forthe first term in Eq. (7) and O(N) for the second term (here, O is theLandau symbol). The sum of both terms are summarized to O(N²). Incontrast, if the method of this embodiment is used, because theconvergence of the solution is high, the energy does not need to becalculated, and the amount O(N²) of computation can be saved.

On the other hand, the method of this embodiment rearranges the spins inStep 111. The amount of the processing is estimated as follows. If eachspin is simply compared with the other spins in a repeated manner, theamount is O(N²). This is the upper limit for the processing. The lowerlimit is estimated as follows. Let us consider that an array of spinshas already put in the descending order, and that only the adjacentspins are compared and exchanged. If there is no exchanging, the amountof computation is O(N) because each spin is compared with only one sideof adjacent spins. Because a huge number of exchanges is generally rare,the actual amount of computation approaches O(N). Thus, the overhead ofthis embodiment is about O(N), which is smaller than the overhead O(N²)of the first embodiment.

The amount of computation for the entire local-field response method isdominated by the computation of the effective magnetic field in Step 102a. Because every N site amounts to O(N), the total amount is O(N²).Therefore, the overhead O(N) in this embodiment is negligible in asystem in which N is sufficiently large.

As described above, in the example of FIG. 2, Steps 102 a and 103 areprocessed independently and in parallel for all sites. For this reason,although the processing speed is high, quantum entanglement is notsufficiently taken into consideration. In contrast, in the example ofFIG. 4 in this embodiment, the effect of quantum entanglement isphenomenologically incorporated through adding Steps 112 and 111, whichmakes the influence of the other spins for the relevant spin at themoment reflect at steps 102 a and 103. As a result, solution accuracy isimproved; the convergence of the solution is improved; and thecomputation time can be saved.

Third Embodiment

Quantum-mechanically, the effective magnetic field is determined basedon Eq. (4). An eigenvalue of σ̂_(k) ^(z) is ±1. However, because thelocal-field response method operates such that a spin variable s_(k)^(z) takes an expectation value <σ̂_(k) ^(z)>, |s_(k) ^(z)|≤1 issatisfied. For this reason, the term Σ_(k(≠j))J_(kj)s_(k) ^(z) isgenerally underestimated compared with g_(j).

If the computation is performed while the term of Σ_(k(≠j))J_(kj)s_(k)^(z) is underestimated, the solution accuracy is degraded. Therefore,the value of g_(j) is normalized with reference to the value of s_(k)^(z). A factor c_(i)=(Σ_(k)s_(k) ^(z)(t_(i−1))²/N)^(1/2) is multipliedto g_(j) to obtain g_(j) ^(norm)(t_(i))=c_(i)g_(j). If g_(j)^(norm)(t_(i)) is set as a local term, the contributions of the termsg_(j) ^(norm)(t_(i)) and Σ_(k(≠j))J_(kj)s_(k) ^(z) are almost equal, andthe solution accuracy is improved. Here, let m (t_(m)≤τ) be the numberof divisions in the discrete time axis, and c₁ is set as about c₁=1/m.If c₁ is simply determined in accordance with c_(i)=(Σ_(k)s_(k)^(z)(t_(i−1))²/N)^(1/2) and s_(k) ^(z)(t₀)=0, then c₁=0. The setting ofc₁=1/m is to cope with c₁=0.

FIG. 5 illustrates a flowchart containing the above processes. Adifference from FIG. 4 is that Step 102 a is replaced with Step 102 b.In Step 102 b, the process about the factor c_(i)=(Σ_(k)s_(k)^(z)(t_(i−1))²/N)^(1/2) is added.

FIG. 6 illustrates a modification. If a parameter c_(a) is introducedand the factor c_(i) is set as c_(a)·c_(i), the magnitude of the factorc_(i) can be adjustable. FIG. 6 is such a case. As described below, thefactor c_(i) in Step 10 c in FIG. 9 is also replaced with c_(a)·c_(i).The value of c_(a) is empirically determined, and is about 1 to 50.

Fourth Embodiment

In quantum-mechanical spin systems, the spins always affect other spinswith each other. That is, spin σ̂_(j) ^(z) at a site j affects spineσ̂_(k) ^(z) at another site k, and inversely σ̂_(k) ^(z) affects σ̂_(j)^(z). Therefore, the spin σ̂_(j) ^(z) affects itself through the spinσ̂_(k) ^(z) at site k. That is the reason why, in quantum mechanics, thestate of a spin depends not only on the state of other spins interactingwith the relevant spin but also on the state of itself. A magnitude ofthe influence on itself through the interaction is proportional toΣ_(k(≠j))J_(kj) ². The first embodiment has described the averagedquantum effect. Because Σ_(k(≠j))J_(kj) ² are squared terms, they areleft even after averaged. That is the reason why the averaged responsefunction becomes r_(b) ⁰(t)≠1. The open circles and the solid line inFIG. 3 show such a behavior. A detailed theory about Σ_(k(≠j))J_(kj) ²is disclosed in NPL 2.

Terms Σ_(k(≠j))J_(kj) ² that are the origin of r_(b) ⁰(t)≠1 containinformation J_(kj) for each problem. If the information is usable, thecomputation becomes more accurate. Thus, we replace the averagedresponse function r_(b) ⁰(t) with r_(b) ⁰ _(mod)(t), where r_(b) ⁰_(mod)(t) is defined by 1−r_(b) ⁰ _(mod)(t)=(1−r_(b)⁰(t))Σ_(k(≠j))J_(kj) ²/Σ_(k(≠j))ave(J_(kj) ²). Herein, ave(J_(kj) ²) isthe average of J_(kj) ² which are used to determine r_(b) ⁰(t). With theimprovement, the response function reflects an actual problem, and thesolution accuracy increases.

FIG. 7 illustrates a flowchart containing the above processes. Adifference from FIG. 5 is that Step 103 is replaced with Step 103 c.

Fifth Embodiment

Quantum entanglement and linear superposition are representativecharacteristic natures in quantum mechanics. The effect of the formerquantum entanglement has been phenomenologically incorporated in up tothe fourth embodiment. The effect of the latter linear superposition isphenomenologically incorporated in NPL 2. In this embodiment, the botheffects are incorporated at the same time. FIG. 8 illustrates such acase.

In the embodiment illustrated in FIG. 8, the difference from FIG. 7 isthat Step 102 b is replaced with Step 102 d.

The characteristic of linear superposition is prominent in a time regionwhere the sign of a spin changes.

Quantum-mechanically, bands anti-cross in the vicinity of that region,and the state is linear superposition there. As a result, s_(j)^(z)(t)≈0, and correspondingly B_(eff,j) ^(z)(t)≈0. In order tophenomenologically incorporate the effect, the effective magnetic fieldsat times t_(i) and t_(i−1) are linearly combined to give the effectivemagnetic field B_(eff,j) ^(z)(t_(i)) at t=t_(i). Specifically, B_(j)^(z0)(t_(i)) defined by Eq. (8) is first determined using the spinvalues s_(j) ^(z)(t_(i−1)) at time t_(i−1).

$\begin{matrix}{{B_{j}^{z\; 0}\left( t_{i} \right)} = {{\sum\limits_{k{({\neq j})}}{J_{kj}{s_{k}^{z}\left( t_{i - 1} \right)}}} + {c_{i}g_{j}}}} & (8)\end{matrix}$

Next, the term at prior time is taken into consideration, and B_(j)^(z)(t_(i)) defined by Eq. (9) is computed.

B _(j) ^(z)(t _(i))=(1−u)B _(j) ^(z0)(t _(i))+uB _(j) ^(z)(t_(i−1))  (9)

Herein, u is appropriately adjusted within 0≤u≤1 to make the solutionaccuracy high, typically u≈0.1. The effective magnetic field, includinga transverse field and its schedule, is given by Eq. (10).

$\begin{matrix}{{B_{{eff},j}\left( t_{i} \right)} = \left( {{\left( {1 - \frac{t_{i}}{\tau}} \right)\gamma},0,{\frac{t_{i}}{\tau}{B_{j}^{z}\left( t_{i} \right)}}} \right)} & (10)\end{matrix}$

When the effective magnetic field is obtained, s_(j) ^(z)(t) isdetermined using the response function r_(b) ⁰ _(mod)(t) in accordancewith Eq. (11).

s _(j) ^(z)(t _(i))/s _(j) ^(x)(t _(i))=r _(b mod) ⁰ ·B _(eff,j) ^(z)(t_(i))/B _(eff,j) ^(x)(t _(i))  (11)

Sixth Embodiment

Up to the fifth embodiment, we have described the embodiments to improvethe solution accuracy by adding various quantum effects, especially byadding the effect of quantum entanglement. However, the correct solutionis not necessarily derived even if the above methods are used. For thisreason, this embodiment will describe an auxiliary means.

Because the local-field response method is a deterministic method, theresult is always the same if the same process is performed with the sameinitial spin values. If the initial values or the process is changed,the result may be different. Hence, we sweep the magnetic field severaltimes while changing the initial values or the process, and we selectthe spin state giving the lowest energy in the process as the finalsolution to further improve the solution accuracy. FIG. 9 illustrates aflowchart for that case.

In FIG. 9, Processes 10 a to 10 d represent those described in FIG. 2and FIGS. 4 to 8. Energy H_(p) ^(q)(t_(m)) (q=1, 2, . . . ) iscalculated in Steps 303 using spin values s_(j) ^(zfd) obtained inProcesses 10 a to 10 d, and the lowest value E_(best) is determined inStep 304. The spin values giving the lowest energy become the finalsolution.

Process 10 a just performs that in FIG. 2 and FIGS. 4 to 8. In Process10 b, the sign of the result s_(j) ^(zfd) obtained in Process 10 a isinverted, and the sign-inverted value is set as the initial value.Strictly, the initial value is s_(j) ^(z)(t₀)=s_(j) ^(zfd)/div using aparameter div. The reason for dividing it with div is to sufficientlyreduce the initial value, and div is set to be about m. The reason forinverting the sign is that the sign-inverted state is located farthestfrom the uninverted state in the spin-configuration space, and that thesign inversion gives an opportunity to escape from a local minimum if aspin state is trapped there. Process 10 c sets the initial spin valuesas 0 similarly to Process 10 a. However, a factor c_(a) is multiplied toa coefficient c_(i) such as c_(a)·c_(i) for the coefficient of the localterm g_(j) (see the third embodiment) so as to change a balance betweenthe interaction term and the local term. The value of c_(a) isempirically determined, and it is about 1 to 50. Process 10 d determinesthe initial spin values using random numbers. Using random numbersexpands possibilities of finding the solution, and the process can besimplified. Process 10 d is repeatedly performed while changing therandom numbers for the initial values. The solution accuracy increasesas repetitions increase.

Seventh Embodiment

This embodiment is described as an algorithm, which may be executed as asoftware on a general computer or on a dedicated hardware. Thelocal-field response method in the embodiments has a feature that arelatively simple computation is repeated. Therefore, it is effectivethat the repeated computation part is established with a dedicatedhardware, and the other parts are achieved with a general-purposedevice.

FIG. 10 illustrates an example of a configuration for the computation inthis embodiment. FIG. 10 is similar to the configuration of a generalcomputing apparatus, but it includes a local-field response computingdevice 600. The local-field response computing device 600 is thededicated part for the computation described in the first to sixthembodiments. The other general computation is performed with a generalcomputation device 502.

The above configuration may be constructed as a single computer.Alternatively, the configuration may be constructed from differentcomputers connected through a network, where arbitrary parts such as amain storage device 501, a general computation device 502, a controldevice 503, an auxiliary storage device 504, an input device 505, and anoutput device 506 are connected to the network.

General computation is performed in the same procedure as in an ordinarycomputing apparatus. The main storage device 501 (storage unit) and thegeneral computation device 502 (computing unit) repeatedly transfer databetween them so as to progress the computation. At that time, thecontrol device 503 works as a control unit. A program executed with thegeneral computation device 502 is stored in the main storage device 501(storage unit). When the main storage device 501 has an insufficientmemory capacity, the auxiliary storage device 504 that is similarly astorage unit is used. The input device 505 is used for inputting data, aprogram, and the like, and the output device 506 is used for outputtingresults. For the input device 505, not only a manual input device suchas a keyboard but also an interface for a network connection may beused. In addition, the interface serves as the output device as well.

As described in the first to sixth embodiments, N spin variables s_(j)^(z)(t) and N effective magnetic field variables B_(eff,j) ^(z)(t) areiteratively determined in the local-field response computation in theembodiments. The local-field response computing device 600 dedicatedlyexecutes this iterative computation.

In the sixth embodiment, we performed the similar Processes 10 a to 10 dand obtained solutions s_(j) ^(zfd) for the respective processes. Theobtained solutions are transferred from the local-field responsecomputing device 600 to the main storage device 501, and the energyH_(p) ^(q)(t_(m)) and E_(best) are computed using the generalcomputation device 502. That is, individual computations not belongingto the iterated computation is performed using the general computationdevice 502, and the dedication rate in the local-field responsecomputing device 600 is increased.

Eighth Embodiment

High-difficulty problems in combinatorial optimization problems belongto NP-hard. In addition, problems classified into “P” and problemsclassified into “NP” both can be transferred to an NP-hard problem.Therefore, if an NP-hard combinatorial optimization problem is solved,almost all the combinatorial optimization problems are solved. A groundstate-searching problems described Eq. (1) include NP-hard problems. Inthis embodiment, we show how to treat an NP-hard problem in accordancewith Eq. (1) by using a maximum cut (MAX-CUT) problem that is arepresentative NP-hard problem as an example.

The maximum cut problem is a problem of graph theory. In the graphtheory, a graph G is constructed from an node set V and a vertex set E,and is written as G=(V, E). An edge e is written as e={i, j} using twonodes. When a graph is defined by taking the direction of edges e intoconsideration, the graph is called a directed graph. When a graph isdefined without taking the direction into consideration, the graph iscalled an undirected graph. For an edge e, a weight is also defined, andit is written as w_(ij) and w_(ji). For an undirected graph,w_(ij)=w_(ji). Let us think of dividing nodes in a weighted undirectedgraph G=(V, E) into two groups. The MAX-CUT problem is to find adivision method to maximize a total sum of weights of cut edges in thedivision problem. Let G₁=(V₁, E₁) and G₂=(V₂, E₂) be the divided twoundirected graphs. Then, the MAX-CUT problem is to maximize Eq. (12).

$\begin{matrix}{{w\left( V_{1} \right)} = {\sum\limits_{{i \in V_{1}},{j \in V_{2}}}w_{ij}}} & (12)\end{matrix}$

If s_(i)=1 is set for edge i∈V₁, and s_(j)=−1 is set for edge j∈V₂, Eq.(13) is derived.

$\begin{matrix}{{w\left( V_{1} \right)} = {{\sum\limits_{{i \in V_{1}},{j \in V_{2}}}w_{ij}} = {{\frac{1}{4}{\sum\limits_{\underset{({i \neq j})}{i,{j \in V}}}{w_{ij}s_{i}s_{j}}}} = {{\frac{1}{4}{\sum\limits_{\underset{({i \neq j})}{i,{j \in V}}}w_{ij}}} - {\frac{1}{2}{\sum\limits_{\underset{({i \neq j})}{i > j}}{w_{ij}s_{i}s_{j}}}}}}}} & (13)\end{matrix}$

Because the first term in the rightest side is a constant once a graph Gis given, the MAX-CUT problem becomes a problem of minimizingΣ_(i>j)w_(ij)s_(i)s_(j). Because the Hamiltonian of the Ising spin glassmodel is given by Eq. (14),

$\begin{matrix}{{\hat{H}}_{p} = {{- {\sum\limits_{i > j}{J_{ij}{\hat{\sigma}}_{i}^{z}{\hat{\sigma}}_{j}^{z}}}} - {\sum\limits_{j}{g_{j}{\hat{\sigma}}_{j}^{z}}}}} & (14)\end{matrix}$

the MAX-CUT problem is equivalent to the ground state-searching problemof Eq. (14) with J_(ij)=−w_(ij) and g_(j)=0.

Ninth Embodiment

In the second embodiment illustrated in FIG. 4, quantum entanglement istaken into consideration all over the variables at N sites correspondingto N spins in the computation. However, as the value of |s_(j) ^(z)|increases, the probability with which a spin is inverted decreases. Inthis situation, spins may be assumed not to invert after a predeterminedtime, i.e., the effect of quantum entanglement may be ignored after thepredetermined time. In other words, the computation at time t_(i) may bechanged such that the condition of Step 113 in FIG. 4 is changed tol=x<N, and that if “Yes” at Step 113, the spins after x are notsubjected to Step 112, and Steps 102 a and 103 are performedindividually and in parallel for all remaining sites. In this way, it ispossible to shorten a process time while securing accuracy to somedegree.

In this embodiment, the same function as that configured with softwarecan be achieved with hardware such as an FPGA (Field Programmable GateArray) and an ASIC (Application Specific Integrated Circuit).

The invention is not limited to the above embodiments, and variousmodifications can be made. For example, some configurations of a certainembodiment may be replaced with the configurations of anotherembodiment, and the configuration of the other embodiment may also beadded to the configuration of a certain embodiment. In addition, part ofthe configurations of each embodiment may be added to or replaced withpart of the configurations of other embodiments, and some of theconfigurations may be omitted.

What is claimed is:
 1. A computing apparatus which includes a computingunit, a storage unit, and a control unit, and performs computation underthe control of the control unit while transferring data between thestorage unit and the computing unit, wherein N variables s_(j) ^(z)(j=1, 2, . . . , N) take a range of −1≤s_(j) ^(z)≤1, and an assignmentis set with coefficients g_(j) indicating local terms and coefficientsJ_(kj) (k, j=1, 2, . . . , N) indicating inter-variable interactions,time is divided into m, and the computing unit discretely performscomputation from t=t₀ (t₀=0) to t_(m) (t_(m)≤τ), variables B_(eff,j)^(z)(t_(i)) and s_(j) ^(z)(t_(i)) at each time t_(i) (i=1, 2, . . . , m)are determined in this order, B_(eff,j) ^(z)(t_(i)) is a function ofs_(k) ^(z) (t_(i−1)) J_(kj), g_(j), and t_(i), s_(j) ^(z)(t_(i)) is afunction of B_(eff,j) ^(z) (t_(i)) and t_(i), and initial values at timet₀ are set as B_(j) ^(z)(t₀)=0 and s_(j) ^(z)(t₀)=0, for determiningB_(eff,j) ^(z)(t_(i)) and s_(j) ^(z)(t_(i)) at time t_(i) (i=1, 2, . . ., m), first, s_(j) ^(z)(t_(i−1)) are put in descending order such that|s_(m1) ^(z)(t_(i−1))|≤s_(m2) ^(z)(t_(i−1))|≤|s_(m3) ^(z)(t_(i−1))|≤ . .. ≤|s_(mN) ^(z)(t_(i−1))|, then, B_(eff,m1) ^(z)(t_(i)) and s_(m1)^(z)(t_(i)) at site m₁ are determined at the first time, and s_(m1)^(z)(t_(i−1)) is set to be s_(m1) ^(z)(t_(i−1))=sgn(s_(m1)^(z)(t_(i)))|s_(m1) ^(z)(t_(i−1))|, next, B_(eff,m2) ^(z)(t_(i)) ands_(m2) ^(z)(t_(i)) at site m₂ are determined and s_(m2) ^(z)(t_(i−1)) isset to be s_(m2) ^(z)(t_(i−1))=sgn(s_(m2) ^(z)(t_(i)))|s_(m2)^(z)(t_(i−1))|, the computation at site m₃ is performed similarly, andthe computation up to site m_(x) (herein, 3≤x≤N) is performed similarlyfor the computation at time t_(i), and the variable s_(j) ^(z)approaches −1 or 1 as the time step progresses from t=t₀ to t=t_(m), anda solution is determined as s_(j) ^(zfd)=−1 if s_(j) ^(z)<0 and as s_(j)^(zfd)=1 if s_(j) ^(z)>0.
 2. The computing apparatus according to claim1, wherein B_(eff,j) ^(z)(t_(i)) is determined by B_(eff,j)^(z)(t_(i))=(t_(i)/τ)·(Σ_(k(≠j))J_(kj)s_(k) ^(z)(t_(i−1))+g_(j)).
 3. Thecomputing apparatus according to claim 1, wherein B_(eff,j)^(x)(t_(i))=γ(1−t_(i)/τ) is set using a constant γ, θ is define by tanθ=B_(eff,j) ^(z)(t_(i))/B_(eff,j) ^(x)(t_(i)), s_(j) ^(z)(t_(i)) isdetermined by s_(j) ^(z)(t_(i))=sin θ, and thus s_(j)^(z)(t_(i))=f(B_(eff,j) ^(z)(t_(i)), t_(i))=sin{arctan(B_(eff,j)^(z)(t_(i))/B_(eff,j) ^(x)(t_(i)))} is obtained using a function f. 4.The computing apparatus according to claim 3, wherein correctionparameters r_(s) and r_(b) are added to the function f, θ is defined bytan θ=r_(b)·B_(eff,j) ^(z)(t_(i))/B_(eff,j) ^(x)(t_(i)), and s_(j)^(z)(t_(i)) is determined by s_(j) ^(z)(t_(i))=r_(s)·sin θ, and thus,the function f becomes f(B_(eff,j) ^(z)(t_(i)),t_(i))=r_(s)·sin{arctan(r_(b)·B_(eff,j) ^(z)(t_(i))/B_(eff,j)^(x)(t_(i)))}.
 5. The computing apparatus according to claim 2, whereinc_(i)=(Σ_(k)(s_(k) ^(z)(t_(i−1)))²/N)^(1/2) and g_(j)^(norm)(t_(i))=c_(i)·g_(j) are set, and B_(eff,j) ^(z)(t_(i)) isdetermined by B_(eff,j)^(z)(t_(i))=(t_(i)/τ)·(Σ_(k(≠j))J_(kj)s_(k)(t_(i−1))+g_(j)^(norm)(t_(i))).
 6. The computing apparatus according to claim 5,wherein B_(eff,j) ^(z)(t_(i)) is determined by B_(eff,j)^(z)(t_(i))=(t_(i)/τ)·(Σ_(k(≠j))J_(kj)s_(k)(t_(i−1))+c_(a)·g_(j)^(norm)(t_(i))) using a parameter c_(a).
 7. The computing apparatusaccording to claim 4, wherein δr_(b)≡1−r_(b) is defined with respect tothe correction parameter r_(b), and δr_(b) is given asδr_(b)(t)∝Σ_(k(≠j))J_(kj) ².
 8. The computing apparatus according toclaim 5, wherein B_(j) ^(z0)(t_(i))=(Σ_(k(≠j))J_(kj)s_(k)^(z)(t_(i−1))+g_(j) ^(norm)(t_(i))) is defined, B_(j)^(z)(t_(i))=(1−u)B_(j) ^(z0)(t_(i))+uB_(j) ^(z)(t_(i−1)) is definedusing a parameter u satisfying 0≤u≤1, and B_(eff,j) ^(z)(t_(i)) isdetermined by B_(eff,j) ^(z)(t_(i))=B_(j) ^(z)(t_(i))·t_(i)/τ.
 9. Thecomputing apparatus according to claim 1, wherein the computation fordetermining s_(j) ^(zfd) described in claim 1 is performed severaltimes, a parameter div is set to a value as large as m, initial valuesat the second and subsequent computations are set as s_(j)^(z)(t₀)=−s_(j) ^(zfd)/div using the solution s_(j) ^(zfd) for the lastcomputation or are set as s_(j) ^(z)(t₀)=1/div or s_(j) ^(z)(t₀)=−1/divusing a random number, H_(p)=−Σ_(k>j)J_(kj)s_(k) ^(zfd)(t_(i))s_(j)^(zfd)−Σ_(j)g_(j)s_(j) ^(zfd) is calculated for each computation, andthe final solution is s_(j) ^(zfd) giving the minimum H_(p) in therepeated computations.
 10. The computing apparatus according to claim 1,wherein after site m_(x), the computation of time t_(i) is performed forall remaining sites independently and in parallel.
 11. A computingmethod which uses a computing apparatus including a computing unit, astorage unit, and a control unit, and performs a computation under thecontrol of the control unit while transferring data between the storageunit and the computing unit, wherein N variables s_(j) ^(z) (j=1, 2, . .. , N) take a range of −1≤s_(j) ^(z)≤1, and an assignment is set withcoefficients g_(j) indicating local terms and coefficients J_(kj) (k,j=1, 2, . . . , N) indicating inter-variable interactions, time isdivided into m, and the computing unit discretely performs computationfrom t=t₀ (t₀=0) to t_(m) (t_(m)≤τ), variables B_(eff,j) ^(z)(t_(i)) ands_(j) ^(z)(t_(i)) at each time t_(i) (i=1, 2, . . . , m) are determinedin this order, B_(eff,j) ^(z)(t_(i)) is a function of s_(k)^(z)(t_(i−1)), J_(kj), g_(j), and t_(i), s_(j) ^(z)(t_(i)) is a functionof B_(eff,j) ^(z)(t_(i)) and t_(i), and initial values at time t₀ areset as B_(j) ^(z)(t₀)=0 and s_(j) ^(z)(t₀)=0, for determining B_(eff,j)^(z)(t_(i)) and s_(j) ^(z)(t_(i)) at time t_(i) (i=1, 2, . . . , m),first, s_(j) ^(z)(t_(i−1)) are put in descending order such that |s_(m1)^(z)(t_(i−1))|≤|s_(m2) ^(z)(t_(i−1))|≤|s_(m3) ^(z)(t_(i−1))|≤ . . .≤|s_(mN) ^(z)(t_(i−1))|, B_(eff,m1) ^(z)(t_(i)) and s_(m1) ^(z)(t_(i))at site m₁ are determined at the first time, and s_(m1) ^(z)(t_(i−1)) isset to be s_(m1) ^(z)(t_(i−1))=sgn(s_(m1) ^(z)(t_(i)))|s_(m1)^(z)(t_(i−1))|, for the following sites, the same computation isperformed up to site m_(x) (herein, 1≤x≤N) for the computation at timet_(i), and variables s_(j) ^(z) approach −1 or 1 as the time stepprogresses from t=t₀ to t=t_(m), and a solution is determined as s_(j)^(zfd)=−1 if s_(j) ^(z)<0 and as s_(j) ^(zfd)=1 if s_(j) ^(z)>0.
 12. Thecomputing method according to claim 11, wherein the same computation isperformed up to site m_(N) after site m₁ for the computation at timet_(i).
 13. The computing method according to claim 11, wherein aftersite m_(x), all remaining sites are processed independently and inparallel to perform the computation at time t_(i).